

function [widthNormalization,totalIntenNormalization,surroundIntenNormalization,maxIntenNormalization] = widefieldTrajectoryNormalization_2led_v1_0(localgcamplevel,traj,fps)

%% set thresholds/values
gcampthresh = 4;
minStretchLength = 1.5;       %in seconds
averagingLength = 4;        %number of frames to average for normalization quantity
vascularDelay = 1.5;        % in seconds
%% compute some reference values
if nargin < 3 || isempty(fps)
    fps = 1/.14;
end
minStretchLength = round(minStretchLength * fps);
vascularDelay = round(vascularDelay*fps);
numsegs = size(localgcamplevel,1);
numframes = size(localgcamplevel,2);
numstims = size(localgcamplevel,3);

% separate normalizations for each measured trajectory field (not including
% vessel centers)
widthNormalization = nan(numsegs,numstims);
totalIntenNormalization = nan(numsegs,numstims);
surroundIntenNormalization = nan(numsegs,numstims);
maxIntenNormalization  = nan(numsegs,numstims);

for n = 1:numstims
    cur = localgcamplevel(:,:,n) < gcampthresh;
    for n2 = 1:size(cur,1)
        temp = cur(n2,:)';
        i=reshape(find(diff([0;temp;0])~=0),2,[]);      %each column is a stretch of consecutive subthreshold values; top row is start frame, bottom row is end frame
        stretchlengths = diff(i);
        i(2,:) = i(2,:)-1;
        idx = find(stretchlengths > minStretchLength,1);        %use first available stretch
        if isempty(idx)
            continue;
        end
        idx = (i(1,idx) + vascularDelay):(i(1,idx) + vascularDelay + averagingLength - 1);
        if max(idx) > numframes
            continue;
        end
        widthNormalization(n2,n) = median(traj(n).vesselWidth(n2,idx),'omitnan');
        totalIntenNormalization(n2,n) = median(traj(n).vesselTotalInten(n2,idx),'omitnan');
        surroundIntenNormalization(n2,n) = median(traj(n).surroundInten(n2,idx),'omitnan');
        maxIntenNormalization(n2,n) = median(traj(n).maxInten(n2,idx),'omitnan');
    end
end
widthNormalization = median(widthNormalization,2,'omitnan');
totalIntenNormalization = median(totalIntenNormalization,2,'omitnan');
surroundIntenNormalization = median(surroundIntenNormalization,2,'omitnan');
maxIntenNormalization = median(maxIntenNormalization,2,'omitnan');
